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Abstract 

In  the  winter  of  1997,  a  shallow  water  experiment  was  conducted  near  the  Sctipps  Pier  in  La 
Jolla,  California,  USA.  The  goal  was  to  determine  the  dynamics,  distribution  and  acoustic  effects  of 
bubbles  just  offshore  from  active  surf.  A  major  component  of  the  experiment  was  the  “Delta 
Frame,”  an  apparatus  that  supported  two  acoustic  sources  and  eight  receivers.  Acoustic  intensity 
was  measured  at  frequencies  between  39  and  244  kHz  over  the  resulting  16  horizontal  ray  paths. 
Paths  ranged  in  length  from  2.5  to  8.6  m.  In  the  present  paper,  a  tomography  algorithm  is 
developed  and  implemented  using  Delta  Frame  data.  Measurements  are  combined  to  produce 
quantitative  cross-sectional  images  of  the  attenuation  associated  with  the  bubble  cloud.  Numerical 
simulations  predict  that  the  frame  can  resolve  details  of  the  field  down  to  about  2  m.  Images 
constructed  at  different  acoustic  frequencies  are  scaled  and  compared.  A  five-minute  sequence  of 
images  is  studied  in  detail.  Swell  waves  are  shown  to  cause  rapid  fluctuations  in  the  images  due  to 
the  vertical  structure  of  the  bubble  cloud.  At  a  wave  trough,  the  bubble  cloud  dips  into  the 
horizontal  plane  probed  by  the  Delta  ftame  and  the  attenuation  increases. 


Introduction 

Bubbles,  even  in  small  concentrations,  have  strong  effects  on  underwater  acoustic 
propagation.  They  produce  strong  attenuation,  and  they  change  the  sound  speed.  Unlike  most 
acoustic  environments,  bubble  clouds  cause  the  sound  to  be  strongly  dispersed.  As  a  crude 
approximation  useful  to  guide  one's  thinking,  the  complex  index  of  refraction  at  any  frequency 
depends  on  the  bubbles  that  resonate  at  that  frequency,  and  the  resonant  frequency  is  inversely 
proportional  to  the  radius  of  the  bubble.  The  processes  that  remove  bubbles  from  the  water 
determine  the  index  of  refraction  as  a  function  of  frequency.  Large  bubbles  rapidly  rise  to  the 
surface,  and  small  bubbles  dissolve,  leaving  intermediate  size  bubbles  to  dominate  the  acoustic 
effects. 


distribution  statement  a 

Approved  for  Public  Release 
Distribution  Unlimited 

me  QUALiry  inspected  a 


20000313  066 


Because  of  the  strong  dispersion,  a  transnnitted  pulse  shape  is  greatly  distorted  if  it 
propagates  a  significant  distance  through  a  bubble  cloud.  Under  such  conditions,  the  very 
definition  of  a  travel  time  becomes  uncertain.  It  turns  out  to  be  much  easier  to  interpret  intensity 
measurements  in  terms  of  attenuation  rather  than  travel  times.  The  real  part  of  the  index  of 
refraction  is  better  studied  by  short-range  phase  velocity  measurements  than  by  transmission 
(group)  speeds.  Another  consequence  of  the  strong  dispersion  is  that  measurements  at  different 
frequencies  give  independent  information  about  the  bubble  clouds. 

Near  the  ocean  shore,  bubbles  produced  by  surf  are  often  the  dominant  complication  in 
acoustic  propagation.  These  bubbles  are  transported  by  nearshore  currents.  Rip  currents  carry 
them  offshore,  and  may  affect  attempts  to  use  acoustics  probes  near  the  shore.  An  experiment  to 
determine  the  dynaniics,  distributions,  and  acoustic  effects  of  bubbles  in  shallow  water  just  offshore 
from  active  surf  was  conducted  from  24  February  to  14  March,  1997.  The  region  of  interest  was 
just  north  of  the  pier  at  the  Sciipps  Institution  of  Oceanography.  The  region  extended  from  a  few 
tens  of  meters  to  a  thousand  meters  from  the  beach.  Efforts  were  made  to  acquire  a  comprehensive 
data  set  of  coastal  dynamics,  bubble  distributions,  and  acoustic  propagation  using  a  variety  of 
techniques  [1-5].  A  major  emphasis  was  placed  on  a  smaller  region  approximately  400  m^  in  area 
and  about  300  m  offshore,  where  most  of  the  instruments  were  clustered. 

The  dynamics  of  rip  currents  at  the  Scripps  Pier  have  been  discussed  in  a  number  of  papers 
[6-8].  An  earlier  sector-scan  sonar  study  of  tips  at  the  SIO  pier  found  clear  evidence  of  tip  currents 
and  described  their  character  which  includes  strong  acoustic  backscatter  [8].  That  study  suggested 
that  the  backscatter  is  due  to  inicrobubbles  rather  than  the  sediment  load  they  also  carry,  but  could 
not  confirm  the  speculation.  The  bubble  content  of  tip  events  was  reported  to  be  one  of  their  least 
understood  characteristics.  A  sector-scan  sonar  used  during  the  Scripps  Pier  Bubble  Experiment 
also  observed  the  acoustic  backscattering  [3].  A  preponderance  of  evidence  collected  by  the 
participants  established  the  fact  that  backscattering  is  due  to  microbubbles  [1-5,9,10]. 

An  important  component  of  thel997  experiment  was  the  Delta  Frame.  The  frame  supported 
two  sources  and  eight  hydrophones.  Intensity  (interpreted  by  attenuation)  and  timing  (to  try  to 
extract  travel  time)  measurements  were  made  along  the  resulting  sixteen  intersecting  paths  at  eight 
frequencies  between  39  and  244  kHz.  This  fiequency  range  spanned  the  resonant  frequencies  of 
bubbles  thought  to  be  present.  Note  that  individual  transmission  measurements  represent 
intrinsically  averaged  quantities.  The  logarithm  of  intensity  (i.e.,  the  signal  level  in  dB),  for 
example,  is  reduced  by  attenuation  integrated  over  the  entire  acoustic  path.  In  order  to  recover 
details  of  the  field  smaller  than  the  path  length,  it  is  necessary  to  combine  measurements  taken 
along  different  paths.  Tomography  is  a  systematic  way  to  combine  data  and  produce  quantitative 
cross-sectional  pictures  of  the  medium. 


In  this  paper,  a  tomography  algorithm  is  first  derived  and  then  applied  to  the  Delta  Frame 
measurements.  Features  of  the  bubble  field  between  2  and  10  m  in  size  are  imaged.  The  images 
show  how  the  bubble  field  responds  to  ensonification  at  different  acoustic  frequencies  and  how  the 
field  evolves  in  time.  Modulation  of  the  bubble  field  by  local  surface  waves  is  observed. 

In  Section  1,  relevant  details  of  the  Delta  Frame  measurements  are  discussed.  Section  2 
begins  with  a  derivation  of  the  tomography  algorithm.  The  algorithm  is  then  implemented  and 
tested  in  numerical  simulations  to  estimate  the  resolution  that  can  be  expected.  The  algorithm  is 
applied  to  the  Delta  Frame  data  in  Section  3.  The  frequency  dependence  of  attenuation  is 
examined.  A  sequence  of  images  sampling  the  bubble  field  every  second  is  presented,  finally,  the 
results  are  discussed  in  Section  4. 

L  Delta  Frame  Measurements 

The  Delta  Frame  was  deployed  approximately  12  m  north  of  the  Scripps  Rer  near  piling 
33.  The  outer  frame  was  in  the  shape  of  a  triangle  with  each  leg  approximately  9.4  m  long.  Figure 
1  is  a  top  view  showing  the  location  of  the  two  acoustic  sources  labeled  A  and  B  at  vertices  of  the 
triangle.  The  eight  numbered  receivers  are  distributed  around  the  outer  frame.  Also  shown  are  the 
resulting  16  intersecting  ray  paths  that  range  in  length  from  2.5  to  8.6  m.  Legs  1.6  m  in  length 
supported  the  frame.  To  avoid  spurious  reflections,  extenders  supported  the  sources  and  receivers 
above  the  frame  boosting  the  acoustic  paths  to  2.8  m  above  the  seabed  in  water  nominally  4.3  m 
deep.  A  weak  ^vailing  current  to  the  south  insured  that  effects  of  the  pier  pilings  did  not 
contaminate  the  bubble  fields  at  the  frame.  The  sources  had  narrow  beampattems  in  the  vertical 
directions  to  further  reduce  the  effects  of  scattering  from  the  frame  and  to  eliminate  scattering  from 
the  surface  and  bottom. 

Measurements  along  the  paths  were  made  at  eight  frequencies:  39, 59, 78, 98,  117,  146,  186, 
and  244  kHz.  Each  frequency  component  had  a  0.1  ms  duration,  with  a  modified  Hanning 
envelope.  The  fiequency  components  and  the  sequence  of  their  emissions  were  specially  coded  to 
insure  that  the  various  components  did  not  overly  at  any  hydrophone  and  that  all  transmissions 
would  occur  within  a  12  ms  period  to  fit  into  an  active  ping  schedule  including  emissions  from 
other  experiments.  The  received  signals  were  sampled  at  625  kHz  in  12  bits.  The  cycle  of  eight 
transmissions  from  each  of  the  two  sources  were  made  once  a  second  for  an  hour  and  a  half,  this 
constituting  a  Run.  Over  the  course  of  ten  days,  ten  mns  were  made. 

The  conditions  that  prevailed  for  most  of  Run  6  provided  a  very  stable  baseline  of  signals. 
For  Source  A  the  standard  deviation  of  the  signal  level  was  less  than  0.02  dB  and  for  Source  B  less 
than  0.05  dB  during  quiescent  periods.  Over  the  course  of  the  ten  days  of  measurements,  numerous 
periods  occurred  that  matched  the  data  of  Run  6  within  0.1  dB  (equating  to  errors  of  less  than  0.05 
dB/m)  and  no  period  occurred  in  which  data  exceeded  the  data  of  Run  6  by  more  than  the  expected 


error.  Consequently,  averaging  ten  pings  in  Run  6  beginning  with  Sequence  0200  establishes  a 
valid  "zero-bubble"  baseline.  This  baseline  calibrates  the  "no  loss"  case  and,  unless  otherwi^  noted, 
all  values  of  loss  are  referenced  to  it 

These  small  variations  in  transmission  properties  are  to  be  compared  to  the  attenuation  due 
to  bubbles  that  were  determined  to  be  5  to  30  dB  when  rip  events  prevailed.  Generally  these  events 
were  more  prevalent  during  low  tide  when  surf  breaking  was  just  shoreward  of  die  frame.  The 
remainder  of  this  work  will  discuss  these  variations  in  attenuation. 

11.  Tomographic  Inversion  Algorithm 

It  is  not  practical  to  image  individual  bubbles.  Rather,  the  goal  of  the  inversion  is  to  determine  the 
effective  index  of  refraction  produced  by  an  aggregate  of  bubbles.  The  effective  index  of  refraction 
is  complex  having  a  real  part  related  to  the  speed  of  sound  and  an  imaginary  part  related  to 
attenuation.  In  the  following,  the  inversion  algorithm  is  derived.  The  algorithm  is  then  implemented 
and  tested  in  numerical  simulations  to  estimate  the  resolution  of  the  imaging  system. 

2.a  Theory 

Let  a(r)  represent  the  spatially  varying  attenuation  associated  with  the  bubble  field  in  a  horizontal 
plane  r  =  (x,y).  Sound  is  transmitted  through  the  bubble  field  from  several  different  directions 

within  the  plane.  Taking  the  negative  of  the  logarithm  of  the  intensity  measured  along  path  i  yields 
the  datum  4  modeled  by 

</,  =Ja(r,)ifn,.  +  e,.  (1) 

The  vector  r,  is  constrained  to  lie  along  the  ray  path,  is  the  increment  along  the  path,  and  e,.  is 
the  noise.  In  the  present  application,  the  ray  path  is  assumed  to  be  straight  connecting  the  source 
and  receiver.  It  is  assumed  that  N  such  measurements  are  available  probing  the  medium  from 
different  directions.  The  problem  reduces  to  estimating  the  unknown  two-dimensional  function 
a(r)  from  N  measured  projection  integrals. 

The  reconstruction  of  a  function  from  its  projections  is  a  classic  problem  with  wide 
application  in  many  fields.  In  a  medical  CT  scaimer,  for  example,  upwards  of  one  millicm 
projections  are  measured  and  used  to  generate  an  image  [11].  In  such  cases,  well-known 
techniques  based  on  the  Radon  transform  can  be  applied.  These  methods  build  essentially  no 
assumptions  about  the  field  into  the  inversion  algorithm;  the  large  number  of  measurements 
obviates  the  need  for  a  priori  information.  In  the  present  application,  N  equals  only  16  and  an 
alternative  inversion  approach  is  required. 


The  recent  monograph  by  Munk  et  al.  offers  a  comprehensive  review  of  ocean  acoustic 
tomography  [12].  Given  a  limited  number  of  measurements,  the  most  common  approach  is  to 
parameterize  the  medium.  This  is  the  so-called  discrete-data/discrete-solution  approach.  The 
medium  is  divided  into  layers  or  cells  and  the  goal  of  the  inversion  is  to  estimate  the  associated 
parameters.  In  the  present  application,  the  medium  is  not  segmented;  rather,  a  discrete- 
data/continuous-solution  is  sought.  This  approach  is  more  commonly  used  in  geophysics.  The 
continuous-solution  method  produces  an  image  that  is  smooth  and  consistent  with  the  assumed 
statistical  nature  of  the  bubble  field. 

An  estimate  c((r)  for  the  attenuation  field  is  sought  in  the  form 

«(■•)  =  Xft(rK,  (2) 

1=1 

where  the  functions  g,(r)  are  to  be  determined.  Equation  (2)  produces  an  estimate  that  is  linearly 
dependent  on  the  data. 

There  are  several  ways  to  derive  essentially  equivalent  inversion  algorithms.  Often  the 
discussion  couched  in  terms  of  the  Gauss-Markov  method  [12].  The  resulting  image  is  interpreted 
as  a  realization  of  a  random  process.  In  the  present  case,  some  of  the  requirements  of  the  Gauss- 
Markov  model,  specifically  that  both  the  field  a(r)  and  the  noise  be  realizations  of  Gaussian 
processes,  need  not  be  satisfied.  Only  when  they  are  Gaussian,  though,  does  (2)  give  the  optimum 
image  in  the  least-squares  sense.  What  is  required  is  that  both  the  field  and  the  noise  have  known 
autocorrelation  functions  and  that  the  processes  are  uncorrelated  with  each  other.  Although  not 
necessary,  for  simplicity  in  the  development  it  is  further  assumed  that  the  field  is  statistically 
homogeneous  and  isotropic.  Using  angle  brackets  to  denote  an  ensemble  average,  the 
autocorrelation  of  the  attenuation  field  is  defined  by 

(a(r)a(r'))  =  i?(|r-r'|).  (3) 

Define  the  objective  function  J  from  the  difference  between  the  true  field  and  its  estimate: 

7  =  ([a(r)-a(r)f).  (4) 

The  functions  g,.(r)  are  found  by  minimizing  (4).  Let  g  and  d  be  column  vectors  with  elements  gi 
and  di,  respectively.  Substituting  (2)  into  (4)  and  differentiating  gives 


^  =  -(2[a(r)-g'-(r)d]d'-). 

(5) 

Equating  to  zero  and  solving  yields 

.g''(r)  =  (a(r)d'-).(dd'-)'‘ . 

(6) 

The  uncorrelated  nature  of  the  field  and  the  noise,  together  with  (1),  permits  the  ensemble  averages 
in  (6)  to  be  calculated: 

ig.W  =  (o(r)d,)= Js(|r-  . 

(7) 

(dd^)  =  M  +  V, 

(8) 

where  M  is  the  "influence  matrix"  describing  the  correlation  between  pairs  of  measurements.  It  has 
elements 

(9) 

and  V  describes  the  noise  correlation  with  elements  V^.  =  )  • 

Using  (6)-(9),  the  inversion  algorithm  can  be  rewritten  in  a  more  useful  form 

a(r)  =  2^,V.(r),  (10) 

! 

where  the  weights 

q  =  (M  +  V)-*d.  (11) 

A  simple  interpretation  can  be  attached  to  the  inversion  algorithm  (lO)-(ll).  In  the 
measurement  (forward)  problem,  data  are  taken  along  discrete  ray  paths.  In^licitly,  each  ray  is  of 
infinitesimal  width.  The  functions  \ir,.(r)  are  deflned  as  the  backprojection  rays  used  in  the 
reconstruction  (inverse)  problem.  From  its  definition  (7),  each  tjr,.(r)  represents  the  corresponding 
original  ray  path  in  the  forward  problem  convolved  with  the  autocorrelation  of  the  medium.  Data 
taken  along  isolated  rays  in  the  forward  problem,  then,  are  backprojected  over  a  smeared  region  in 
the  inverse  problem.  TTie  precise  nature  of  this  smearing  depends  upon  the  assumed  correlation 
function  of  the  underlying  medium  that  is  being  reconstructed.  Ibe  r-th  backprojection  ray  is 
weighted  by  the  associated  qi .  The  weighted  backprojection  rays  are  then  superimposed  via  (10)  to 


yield  a  continuous  reconstruction.  The  result  is  analogous  to  the  well  known  filtered  backprojection 
method  used  when  there  is  continuous  data  [1 1]. 

Consider  the  computational  burden  associated  with  implementing  the  inversion  algorithm. 
Elements  inthe  NxN  symmetric  influence  matrix  M  are  found  by  calculating  two-dimensional 
integrals.  For  realistic  configurations,  these  integrals  will  likely  need  to  be  calculated  numerically. 
Note,  however,  that  M  is  independent  of  the  measurements;  it,  as  well  as  the  inverse  of  (M+ V)  , 
need  be  calculated  only  once  even  if  new  data  d  are  measured.  Similarly,  each  \(t,  (r)  is  independent 

of  the  data.  Consequently,  all  of  the  computationally  expensive  work  can  be  completed  prior  to 
making  any  measurements.  As  new  measurements  are  made,  either  at  a  later  time  or  at  a  different 
acoustic  frequency,  new  coefficients  are  found  by  the  simple  matrix  multiplication  in  (11).  The  new 
coefficients  are  applied  to  the  precalculated  backprojection  rays  \|r,(r)  and  combined  via  (10)  to 

yield  the  new  image.  Although  it  was  not  implemented  during  the  Scripps  Pier  experiment,  there  is 
clearly  the  potential  to  do  nearly  real  time  imaging. 

2.b  Implementation 

To  evaluate  the  integrals  in  (7)  and  (9),  it  is  convenient  to  define  a  rotated  coordinate  system: 

t  =  xcosG;  +  ysin  0, 

.  (12) 

T1  = -xsinO, -i-ycos0, 

In  the  (^,11)  coordinates,  the  ray  path  lies  parallel  to  the  q-axis,  displaced  by  the  distance 
(Figure  2).  The  source  is  at  q,,  the  receiver  at  q^.  Equation  (7)  becomes 

Vi(*')  =  Jl?([(4-U^+(q-q,-)^F)rfq,  .  (13) 

ii. 

It  has  been  assumed  in  the  derivation  that  the  autocorrelation  of  the  field  is  known.  Because 
the  tme  autocorrelation  function  is  not  known  in  practice,  a  functional  form  must  be  selected  largely 
on  the  basis  of  ^  a  priori  guess  as  to  the  shape  of  the  bubble  cloud.  After  the  reconstraction  is 
completed,  the  parameters  of  the  assumed  function  can  be  varied  to  confirm  that  the  image  is  not 
unduly  sensitive  to  the  initial  choice.  It  is  also  desirable  that  the  candidate  correlation  function  has 
convenient  mathematical  properties.  Accordingly,  a  Gaussian  correlation  function  is  selected: 


/?(|r  -  r'l)  =  exp(-|r-  rf/Lj), 


(14) 


where  Lg  is  the  correlation  length  and  is  the  single  free  parameter  of  the  model.  Substituting  (14) 
into  (13)  yields 


where  the  remaining  integral 


Jexp[-(ti  - 

=  |{erf[(»l,  -  r|)/i,)  +  erf((n  -  ti,)/!,]}. 


(16) 


is  evaluated  in  terms  of  error  functions.  Note  that  the  backprojection  ray  is  centered  about  the 
original  at  ^  =  ^,.  The  width  is  determined  by  L^.  Along  the  direction  of  propagation,  the 

properties  are  dictated  by  (16).  Let  L  =  ti,  -  -p,  be  the  length  of  the  original  ray.  Figure  3  shows 
W(ti;Z,/2,-I/2)  for  various  ratios  of  the  path  length  to  the  correlation  length.  When  L»  L^, 

W  ~l  over  most  of  the  path,  falling  to  1/2  at  the  source  and  receiver  locations,  and  then  rapidly  to 
zero.  As  the  ratio  is  decreased,  the  backprojection  ray  becomes  more  smeared. 

A  final  detail  on  the  implementation:  from  Fig.  1,  it  is  apparent  that  the  difference  between 
measurements  oil  path  A2  and  A1  should  equal  tihe  difference  between  measurements  on  path  B1 
and  B2  to  within  the  experimental  error.  This  redundancy  was  used  to  preprocess  the  data  and 
reduce  the  dimensionality  of  the  inversion. 


2.C.  Numerical  Simulations 

To  better  understand  the  properties  of  the  reconstruction  algorithm  before  applying  it  to  the 
experimental  data,  a  series  of  numerical  simulations  were  performed.  In  this  section,  results  from 
these  simulations  are  summarized. 

A  realization  of  a  two-dimensional  random  process  can  be  synthesized  by  first  generating 
white  noise  and  then  filtering  according  to  the  desired  spectrum.  The  technique  is  conveniently 

implemented  using  fast  Fourier  transforms.  Figure  4(a)  shows  a  10  m  by  10  m  patch  consistent 
with  a  Gaussian  spectrum.  The  correlation  length  is  4  =  2.5  m.  Superimposed  on  the  realization 

is  the  outline  of  the  Delta  Frame.  Not  shown  are  the  N=16  ray  paths  connecting  the  two  sources 
and  eight  receivers  supported  by  the  frame.  Using  two-dimensional  interpolation  and  a  simple 
quadrature  scheme,  the  projection  integrals  along  the  16  paths  can  be  calculated.  This  serves  as  the 
data  that  can  be  input  to  the  inversion  algorithm.  Rgure  4(b)  shows  the  resulting  reconstruction. 
Over  the  interior  of  the  Delta  Frame,  the  results  are  quite  satisfactory  as  all  the  essential  features  of 


the  original  are  recovered  in  the  reconstruction.  Interestingly,  there  is  also  reasonable  agreement 
around  the  periphery  of  the  frame;  it  is  possible  to  extrapolate  to  a  distance  nearly  outside  of  the 

frame.  This  is  due  to  where  the  sources  and  receivers  are  positioned  on  the  frame  and  because  the 
finite  correlation  length  of  the  medium  is  known.  As  one  would  expect,  no  information  about  the 
field  far  outside  of  the  frame  can  be  gleaned. 

Simulations  of  the  type  shown  in  Rg.  4  were  repeated  for  various  values  of  L^.  To  quantify 

the  inversion  quality,  the  difference  between  the  original  and  the  reconstruction  was  calculated  on  a 
pixel-by-pixel  basis.  This  error  was  then  squared  and  averaged  over  the  interior  of  the  delta  frame. 
The  root  mean  square  error  was  normalized  so  that  if  the  original  and  the  reconstruction  were 
uncprrelated  the  error  would  be  one.  Perfect  agreement  would  result  in  zero  error.  Figure  5  «hows 
the  normalized  rms  error  plotted  Versus  the  correlation  length.  As  one  would  expect,  when  the 
correlation  length  is  large  the  error  is  small.  When  the  correlation  length  becomes  less  than  about  2 
m,  the  error  increases  rapidly.  Physically,  this  means  that  there  are  not  enough  ray  paths  to  recover 
the  detailed  structure  in  a  medium  with  scale  sizes  less  than  a  few  meters. 

in.  Reconstruction  of  the  Bubble  Field 

In  this  section,  the  tomographic  inversion  algorithm  is  applied  to  attenuation  measurements  made 
during  the  Scripps  Pier  experiment.  All  the  images  are  oriented  to  be  consistent  with  the  geometry 
shown  in  Fig.  1.  The  effect  of  varying  the  assumed  correlation  length,  the  lone  free  parameter  of 
the  inversion  algorithm,  is  studied.  The  frequency  dependence  of  attenuation  is  demonstrated.  An 
extended  sequence  of  images  is  generated  showing  the  temporal  evolution  of  the  bubble  field. 

Recall  that  each  transmission  sequence  consisted  of  eight  frequencies  and  was  repeated 
every  second.  As  a  labeling  convention,  sequence  number  is  equivalent  to  time  in  seconds  since  the 
start  of  a  run.  There  were  ten  runs  at  Scripps  Pier  Experiment,  each  lasting  approximately  90 
minutes.  Results  in  this  section  concentrate  on  Run  5  which  began  7  March,  1997  at  13:58:20  local 
time,  and  Run  7  which  began  14:23:32  on  8  March.  Run  5  is  of  particular  interest  since  this  is  both 
a  period  of  high  activity  and  a  time  when  other  investigators  were  making  independent 
measurements  of  the  bubble  field. 

3.a  Effect  of  varying  the  correlation  length 

Figure  6(a)  shows  a  sample  bubble  field  reconstruction  from  Run  7.  The  59  kHz  measurements 
were  averaged  over  ten  seconds  beginning  2823  seconds  into  Run  7.  A  correlation  length  ^  =  2.5 

m  was  assumed.  The  spatially  varying  nature  of  the  medium  is  apparent;  a  swath  of  relatively  high 
attenuation,  7  dB/m,  is  visible  over  the  interior  of  the  finme. 

As  formulated,  the  single  free  paranieter  in  the  inversion  algorithm  is  the  correlation 
length  Reducing  the  assuming  correlation  length  reduces  the  overlap  between  pairs  of 


backprojection  rays  (7)  and  forces  the  off-diagonal  terms  in  the  influence  matrix  (9)  to  shrink.  The 
value  for  Lg  was  varied  to  test  the  effect  on  the  inversion.  Figure  6(b)  is  generated  using  the  same 

data  as  Fig.  6(a)  but  assumes  =  1,0  rti.  Reducing  the  correlation  length  has  made  the  image 

more  fragmented  as  it  becomes  concentrated  in  the  areas  where  the  rays  passed  in  the  forward 
problem.  Still,  the  difference  between  the  two  images  is  not  too  great;  regions  of  high  attenuation 
are  located  at  the  same  place  and  are  of  essentially  the  same  level. 

There  are  trade-offs  in  selecting  the  assumed  correlation  length.  Choosing  a  short 
correlation  length  reduces  the  residual  error  between  the  measurement  and  the  least-squares  fit,  but 
at  the  expense  of  a  fragmented  image  like  Fig,  6(b).  Choosing  a  correlation  length  that  is  too  large, 
.however,  can  introduce  oscillations  in  the  image.  This  is  illustrated  in  Fig;  7  using- -data  from 

sequence  3159  of  Run  5.  Figure  7(a)  shows  the  reconstruction  using  the  146  kHz  measurements 
and  =2.5  m.  A  region  where  the  attenuation  takes  a  significantly  negative  value  is  apparent. 

This  non-physical  result  can  be  corrected  by  reducing  the  correlation  length.  In  Fig.  7(b),  =  0.5 

m  was  assumed.  Now  a  small  region  of  high  attenuation  is  seen  occupying  the  space  between 
Source  A  and  Receiver  3.  There  is  essentially  no  other  activity  in  the  image  space  and  the 
reconstruction  goes  only  slightly  negative.  Clearly,  this  is  a  situation  where  using  our  formulas 
with  a  long  correlation  length  is  not  justified. 

Obtaining  negative  values  for  the  reconstructed  attenuation  serves  as  a  useful  diagnostic  in 
assessing  the  plausibility  of  an  image.  Slightly  negative  regions,  as  in  Fig,  6,  are  accepted  as  an 
intrinsic  feature  of  the  inversion  procedure  and  not  corrected.  For  aberrant  situations,  as  in  Fig. 

7(a),  the  correlation  length  is  reduced  and  a  new  image  generated.  Unless  otherwise  noted, 
=  2.5  m  is  assumed. 

3.b.  Effect  of  varying  the  acoustic  frequency 

The  predominant  mechanism  for  attenuation  in  bubble  fields  is  resonant  scattering.  With  the  Delta 
Frame,  each  sequence  consisted  of  eight  frequencies  allowing  bubbles  of  different  sizes  to  be 
probed  synoptically.  Consequently,  for  each  sequence,  up  to  eight  tomographic  images  can  be 
generated.  These  images  show  how  attenuation  varies  with  acoustic  frequency,  but  also  how 
different  sized  bubbles  are  distributed  spatially. 

Figure  8  shows  bubble  field  attenuation  images  for  acoustic  frequencies  from  39  to  244 
kHz.  The  projection  data  have  been  averaged  over  ten  seconds  starting  from  sequence  2823  of  Run 
7.  Comparing  the  images,  the  most  attenuation  occurs  at  59  kHz.  Using  the  orientation  defined  in 
Fig.  1,  the  59  kHz  image  shows  a  localized  region  of  high  attenuation  reaching  7  dB/m  located 
between  receivers  1  and  2.  A  broad  swath  of  comparable  level  stretches  over  the  body  of  the  frame. 
At  higher  frequencies,  the  basic  shape  of  the  bubble  field  does  not  appear  to  change  significantly, 
but  the  level  of  attenuation  is  reduced. 


Figure  8  demonstrates  the  strong  frequency  dependence  of  attenuation.  It  is  difficult, 
however,  to  compare  the  images;  a  color  bar  appropriate  for  59  kHz  has  too  much  dynamic  range 
for  244  kHz.  To  better  compare  the  results,  it  would  be  useful  to  have  a  model  for  the  functional 
dependence  of  attenuation  on  frequency.  From  this  model,  the  images  could  be  rescaled  and 
compared  using  a  common  color  bar.  The  frequency  dependence  of  attenuation  remains  a  subject 
of  active  research.  In  the  present  context,  an  empirical  scaling  can  be  derived  from  the  images  in 
Fig.  8.  For  a  particular  image  d(r),  define  the  frame-averaged  attenuation 


=  J  6:(r)dr 


(17) 


where  the  integration  is  over  the  triangle  bounded  by  the  Delta  Frame.  Figure  9  shows  d 
calculated  for  each  frame  in  Fig.  8  and  plotted  versus  frequency.  The  maximum  is  at  59  kHz. 
Define  the  frequency-scaled  image  at  a  particular  frequency: 


{scaled  image)  = 


^d(59kHz)^ 

.  «(/)  ; 


x{old  image). 


The  proposed  scaling  was  applied  to  the  images  in  Fig.  8  with  the  results  shown  in  Fig.  10. 

Several  observations  can  be  made  from  Fig.  10.  First,  the  scaling  is  reasonable  up  to  186 
kHz;  the  first  seven  images  essentially  "fit"  on  the  same  color  scale.  A  more  telling  point  is  that  the 
first  seven  images  shown  look  remarkably  similar  to  one  another.  Under  the  assumption  of 
resonant  scattering,  recall  that  large  bubbles  are  sensitive  to  low  frequencies  while  small  bubbles 
respond  to  high  frequencies.  The  similarity  of  the  images  implies  that  in  this  instance  the  bubble 
distribution  is  nearly  homogeneous;  the  small  bubbles  responding  to  186  kHz  sound  are  not 
located  in  a  different  region  than  the  large  bubbles  resonating  at  39  kHz. 

The  proposed  scaling  does,  however,  begin  to  fail  at  the  highest  frequency.  The  244  kHz 
image  actually  has  a  maximum  of  16  dB/m  and  a  minimum  of  -4  dB/m  causing  it  to  saturate  at  both 
ends  of  the  color  scale.  Note,  however,  that  the  location  of  the  maxima  is  generally  consistent  with 
what  was  observed  in  all  the  other  images.  As  in  the  previous  section,  decreasing  the  correlation 
length  would  reduce  the  overshoot  in  the  image.  It  is  unclear  if  the  anomalous  behavior  at  244  kHz 
is  a  genuine  effect  or  an  artifact.  The  relatively  low  attenuation  at  the  highest  frequency  suggests  a 
possible  problem  separating  the  signal  from  noise,  or  errors  in  defining  the  background  attenuation 
level.  Additional  comments  on  the  frequency  dependence  of  the  reconstructions  are  presented  in 
the  following  section. 

3.C  Reconstruction  of  an  extended  time-series 


A  five-minute  portion  from  Run  5  beginning  with  sequence  3050  was  chosen  for  detailed  analysis. 
Figure  11  shows  the  146  kHz  attenuation  field  averaged  over  the  frame.  For  the  first  pairt  of  the 
time  series,  the  average  attenuation  is  less  than  1  dB/m.  Beginning  around  3180  seconds,  the 
average  attenuation  generally  increases  to  sequence  3250,  then  falls  over  the  next  40  seconds.  A 
more  rapid  oscillation  is  superimposed  on  these  broad  characteristics.  The  attenuation  is  negligible 
from  3300  to  3349  s.  this  time  period  was  used  to  correct  for  biases  in  the  individual  paths  and  to 
estimate  the  noise  correlation  matrix  V  in  (8). 

Also  shown  in  Fig.  1 1  is  the  average  attenuation  at  244  kHz.  As  expected,  the  average 
attenuation  is  less  than  was  observed  at  146  kHz.  Note  that  the  two  time  series  generally  track  one 
another. 

During  the  period  of  maximum  activity,  some  paths  were  so  strongly  attenuated  ias  to  be 
unmeasureable  for  frequencies  less  than  146  kHz.  Consequently,  the  imaging  concentrated  on  the 
146  kHz  measurements.  With  a  suitable  model  for  the  frequency  dependence,  these  images  could 
presumably  be  scaled  to  predict  attenuation  characteristics  at  other  fiequencies. 

Figure  12  shows  a  ten  second  sequence  of  images  starting  at  3250  s.  No  attempt  was  made 
to  do  any  time  averaging  of  the  data;  each  frame  was  reconstructed  in  isolation.  Certain  features 
appear  persistent  across  the  images.  Referring  to  Fig.  1  for  the  location  of  the  sources  and 
receivers,  the  broad  feature  extending  from  source  A  to  receiver  5  is  apparent.  The  peak  level 
associated  with  the  feature,  however,  varies  from  frame-to-frame. 

A  possible  cause  of  the  variability  in  the  images  is  wave  motion.  To  test  this  conjecture.  Fig. 
13  shows  a(r),  the  attenuation  averaged  over  the  frame  for  the  two  minutes  of  maximum  activity. 
Superimposed  on  the  plot  is  the  local  water  depth  as  measured  by  a  pressure  sensor  located  near 
source  A  [13].  Comparing  the  two  plots,  attenuation  maxima  correspond  to  wave  troughs  while 
attenuation  minima  correspond  to  wave  crests. 

To  quantify  the  apparent  relationship  between  attenuation  and  surface  height,  a  correlation 
analysis  was  performed.  Figure  14  shows  the  average  attenuation  a(/)  and  a  low-pass  filtered 
vereion  a^(r)of  the  time-series.  The  proposed  model  is  that  the  observed  a(r)  is  a  slowly  varying 
function  modulated  by  the  surface  height  h{t).  To  test  the  model,  the  linear  correlation  coefficient  r 
between  the  ratio  a(/)/a^(r)  and  hif)  was  calculated.  Setting  a(r)>l  dB/m  as  a  threshold  for 
model  validity,  r=-0.75  resulted.  For  a  2  dB/m  threshold,  r=-0.86.  This  strong  anti-correlation 
confirms  that  for  this  section  of  data  the  temporal  variability  is  dominated  by  vertical  motion  of  the 
swell  waves.  The  Delta  Frame  measures  a  two-dimensional,  mostly  horizontal,  stracture.  As  a 
result,  the  vertical  motion  cannot  be  followed  directly;  we  just  see  bubbles  disappear  from  the  plane 
as  the  crest  of  the  wave  approaches,  and  return  as  the  trough  approaches.  During  this  period  of 
high  attenuation,  horizontal  advection  of  the  bubble  field  is  expected,  but  it  is  masked  by  these  wave 
effects. 


IV.  Summary 

The  Delta  Frame  was  designed  to  measure  spatial  and  temporal  structure  of  bubble  clouds  in  a  rip 
current.  When  this  experiment  was  being  planned,  there  was  little  information  on  the  spatial  scales 
of  variability  in  the  bubble  clouds.  The  experiment  was  designed  to  measure  variability  on  scales 
between  2  m  and  10  m.  This  was  a  fortunate  choice,  because  the  tomographic  reconstructions  show 
structure  throughout  this  range.  Our  reconstmcted  images  show  that  there  is  significant  variation 
over  2  m;  possibly  there  are  actual  scales  smaller  than  the  resolution  of  the  delta  frame.  On 
occasion,  it  appears  that  structure  smaller  than  the  resolution  occurs;  setting  the. correlation  length  to 
Im  removes  the  large  negative  bubble  density  from  the  image.  On  the  other  hand,  there  are 
stmctures  in  our  images  that  extend  coiisiderably  more  than  2  m  in  one  direction. 

The  delta  frame  measures  the  bubble  field  in  one  horizontal  plane.  There  is  major  vertical 
structure  that  is  made  apparent  by  the  vertical  advection  from  the  sweD  waves.  This  advective  effect 
dominates  the  temporal  variability  in  the  Delta  Frame  plane.  At  a  wave  trough,  the  bubble  cloud  dips 
into  the  horizontal  plane  probed  by  the  Delta  Frame  and  the  attenuation  increases.  At  a  wave  crest, 
the  cloud  rises  and  attenuation  decreases.  These  wave  effects  tend  to  mask  the  horizontal  advection 
of  the  bubble  cloud. 

TTie  discrete-data/continuous-solution  tomography  algorithm  we  chose  proves  useful.  The 
reconstructed  images  do  not  have  discrete  cells  in  which  the  field  is  assumed  constant;  rather,  the 
structure  of  the  independent  functions  composing  the  images  are  set  by  the  geometry  of  the 
experiment  and  by  the  assumed  correlation  scale.  In  the  algorithm  we  use,  it  is  possible  to  carry  out 
sensitivity  studies  and  to  change  parameters  such  as  the  correlation  length. 
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Figure  captions: 

Figure  1.  Scale  drawing  of  Delta  Frame  showing  location  of  sources  A  and  B  and  eight  receivers. 

Figure  2.  Rotated  coordinate  for  ray  i.  Source  and  receiver  at  and  (^,  ,11;.),  respectively. 

Figures.  Smearing  function  defined  by  (16).  Evaluated  at  =  L/2  and  Ti,  =  -L/2 

where  L  is  the  path  length.  Plotted  for  various  ratios  of  the  path  length  to  the  outer  scale. 

Figure  4.  Reconstruction  simulation,  (a)  Original  image.  Realization  of  process  with  Gaussian 
correlation  function  and  correlation  length  2.5  m.  Triangle  shows  location  of  delta  frame,  (b) 
Tomographic  reconstmction.  Red-green-blue  color  scale  is  arbitrary. 

Figure  5.  Estimating  the  resolution  of  the  imaging  system.  Normalized  error  calculated  by 
comparing  original  and  tomographic  reconstruction.  Results  plotted  versus  correlation  length. 
Error  becomes  large  as  correlation  length  shrinks  below  a  few  meters. 

Figure  6.  Effect  of  choosing  a  correlation  length  on  the  reconstmction  using  experimental  data. 
Run  7,  sequence  2823, 59  kHz  data  with  (a)  Lg=2.5  m,  and  (b)  m. 

Figure  7.  Effect  of  choosing  a  correlation  length  on  the  reconstmction  using  experimental  data. 
Run  5,  sequence  3159, 146  kHz  data  with  (a)  L=2.5  m,  and  (b)  L=Q.5  m. 

Figure  8.  Reconstmction  of  bubble  field  attenuation  at  acoustic  frequencies  from  39  to  244  kHz. 
Run  7,  sequence  2823. 

Figure  9.  Average  attenuation  versus  frequency  for  sequence  2823  of  Run  7.  Calculated  applying 
Eq.  (17)  to  images  in  Fig.  8. 

Figure  10.  Scaled  reconstruction  of  bubble  field  attenuation  at  acoustic  frequencies  from  39  to  146 
kHz.  Run  7,  sequence  2823. 

Figure  11.  Five  minute  sequence  of  attenuation  images  averaged  over  area  of  Delta  Frame.  Both 
146  and  244  kHz  data  used.  Time  in  seconds  relative  to  start  of  Run  5. 

Figure  12:  Ten  second  s^uence  of  reconstmctions.  Time  in  seconds  after  start  of  Run  5.  146 
kl^  data  used  to  generate  images. 

Figure  13:  Comparison  of  reconstracted  attenuation  averaged  over  fi'ame  and  measurement  of  local 
water  depth.  Peaks  in  attenuation  correspond  to  troughs  in  surface  waves.  Averaged  attenuation 
derived  froml46'kHz  data  and  time  relative  to  start  of  Run  5. 

Figure  14.  Attenuation  images  averaged  over  frame  compared  to  lowpass  version  of  data.  Ratio  of 
curves  used  in  text  to  constmct  model  for  effects  of  swell  waves.  Time  relative  to  start  of  Run  5. 


Figure  1.  Scale  drawing  of  Delta  Frame  showing  location  of  sources  A  and  B  and  eight  receivers. 
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Figure  2.  Rotated  coordinate  for  ray  i.  Source  and  receiver  at  (4,,  ti,)  and  respectively. 
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Figure  3.  Smearing  function  W(t];ti^,T|j)  defined  by  (17).  Evaluated  at  T|^  =  L/2  and  Ti,  =  -L/2 
where  L  is  the  path  length.  Plotted  for  various  ratios  of  the  path  length  to  the  outer  scale. 


Figure  4.  Reconstruction  simulation,  (a)  Original  image.  Realization  of  process  with  Gaussian 
correlation  function  and  correlation  length  2.5  m.  Triangle  shows  location  of  delta  frame,  (b) 
Tomographic  reconstruction.  Red-green-blue  color  scale  is  arbitrary. 


Figure  5.  Estimating  the  resolution  of  the  imaging  system.  Normalized  error  calculated  by 
comparing  original  Md  tomographic  reconstruction.  Results  plotted  versus  correlation  length. 
Error  becomes  large  as  correlation  length  shrinks  below  a  few  meters. 


Figure  6.  Effect  of  choosing  a  correlation  length  on  the  reconstruction  using  experimental  data. 
Run  7,  sequence  2823, 59  kHz  data  with  (a)  L^^l.5  m,  and  (b)  ^,=1-0  m. 


Figure  7.  Effect  of  choosing  a  correlation  length  on  the  reconstruction  using  experimental  data. 
Run  5,  sequence  3159, 146  kHz  data  with  (a)  L^=2.5  m,  and  (b)  f-j=0.5  m. 


Figure  8.  Reconstruction  of  bubble  field  attenuation  at  acoustic  frequencies  from  39  to  244  kHz. 
Run  7,  sequence  2823. 
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Figure  9.  Average  attenuation  versus  frequency  for  sequence  2823  of  Run  7.  Calculated  by 
applying  Eq.  (17)  to  images  in  Fig.  8. 
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Figure  10.  Scaled  reconstruction  of  bubble  field  attenuation  at  acoustic  frequencies  from  39  to  146 
kHz.  Run  7,  sequence  2823. 


Figure  11.  Five  minute  sequence  of  attenuation  images  averaged  over  area  of  Delta  Frame.  Both 
146  and  244  kHz  data  used.  Time  in  seconds  relative  to  start  of  Run  5. 


Figure  12:  Ten  second  sequence  of  reconstructions.  Time  in  seconds  after  start  of  Run  5.  146 
kHz  data  used  to  generate  images. 
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Figure  13;  Comparison  of  reconstructed  attenuation  averaged  over  frame  and  measurement  of  local 
water  depth.  Peaks  in  attenuation  correspond  to  troughs  in  surface  waves.  Averaged  attenuation 
derived  froml46  kHz  data  and  time  relative  to  start  of  Run  5. 
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Figure  14.  Attenuation  images  averaged  over  frame  compared  to  lowpass  version  of  data.  Ratio  of 
curves  used  in  text  to  construct  model  for  effects  of  swell  waves.  Time  relative  to  start  of  Run  5. 


